bkmr
Katrina Devick
8/10/2018
Data Manipulation
## read in data and only consider complete data
## this drops 327 individuals, but BKMR does not handle missing data
nhanes = na.omit(read_csv("./Data/studypop.csv"))## Parsed with column specification:
## cols(
## .default = col_double()
## )## See spec(...) for full column specifications.## center/scale continous covariates and create indicators for categorical covariates
nhanes$age_z = scale(nhanes$age_cent) ## center and scale age
nhanes$agez_sq = nhanes$age_z^2 ## square this age variable
nhanes$bmicat2 = as.numeric(nhanes$bmi_cat3 == 2) ## 25 <= BMI < 30
nhanes$bmicat3 = as.numeric(nhanes$bmi_cat3 == 3) ## BMI >= 30 (BMI < 25 is the reference)
nhanes$educat1 = as.numeric(nhanes$edu_cat == 1) ## no high school diploma
nhanes$educat3 = as.numeric(nhanes$edu_cat == 3) ## some college or AA degree
nhanes$educat4 = as.numeric(nhanes$edu_cat == 4) ## college grad or above (reference is high schol grad/GED or equivalent)
nhanes$otherhispanic = as.numeric(nhanes$race_cat == 1) ## other Hispanic or other race - including multi-racial
nhanes$mexamerican = as.numeric(nhanes$race_cat == 2) ## Mexican American
nhanes$black = as.numeric(nhanes$race_cat == 3) ## non-Hispanic Black (non-Hispanic White as reference group)
nhanes$wbcc_z = scale(nhanes$LBXWBCSI)
nhanes$lymphocytes_z = scale(nhanes$LBXLYPCT)
nhanes$monocytes_z = scale(nhanes$LBXMOPCT)
nhanes$neutrophils_z = scale(nhanes$LBXNEPCT)
nhanes$eosinophils_z = scale(nhanes$LBXEOPCT)
nhanes$basophils_z = scale(nhanes$LBXBAPCT)
nhanes$lncotinine_z = scale(nhanes$ln_lbxcot) ## to access smoking status, scaled ln cotinine levels
## our y variable - ln transformed and scaled mean telomere length
lnLTL_z = scale(log(nhanes$TELOMEAN))
## our Z matrix
mixture = with(nhanes, cbind(LBX074LA, LBX099LA, LBX118LA, LBX138LA, LBX153LA, LBX170LA, LBX180LA,
LBX187LA, LBX194LA, LBXHXCLA, LBXPCBLA, LBXD03LA, LBXD05LA, LBXD07LA,
LBXF03LA, LBXF04LA, LBXF05LA, LBXF08LA))
lnmixture = apply(mixture, 2, log)
lnmixture_z = scale(lnmixture)
colnames(lnmixture_z) = c(paste0("PCB",c(74, 99, 118, 138, 153, 170, 180, 187, 194, 169, 126)),
paste0("Dioxin",1:3), paste0("Furan",1:4))
## our X matrix
covariates = with(nhanes, cbind(age_z, agez_sq, male, bmicat2, bmicat3, educat1, educat3, educat4,
otherhispanic, mexamerican, black, wbcc_z, lymphocytes_z, monocytes_z,
neutrophils_z, eosinophils_z, basophils_z, lncotinine_z))
### create knots matrix for Gaussian predictive process (to speed up BKMR with large datasets)
set.seed(10)
knots100 <- fields::cover.design(lnmixture_z, nd = 100)$design
#save(knots100, knots100.PCB, file = "./Supervised/BKMR/saved_model/NHANES_knots100.RData") Fit Models
Group VS Fit with all Exposures using GPP and 100 Knots
### Group VS fit with all exposures using GPP and 100 knots
load("./Supervised/BKMR/saved_model/NHANES_knots100.RData")
##### fit BKMR models WITH Gaussian predictive process using 100 knots
### Group VS fit with all exposures using GPP and 100 knots
set.seed(1000)
## Note: the commented out statement is what generated the saved model fits.
## For this lab we are only going to generate 100 MCMC samples to get a sense
## for what the program outputs.
#fit_gvs_knots100 <- kmbayes(y=lnLTL_z, Z=lnmixture_z, X=covariates, iter=100000, verbose=TRUE, varsel=TRUE,
# groups=c(rep(1,times=2), 2, rep(1,times=6), rep(3,times=2),rep(2,times=7)), knots=knots100)
temp <- kmbayes(y=lnLTL_z, Z=lnmixture_z, X=covariates, iter=100, verbose=TRUE, varsel=TRUE,
groups=c(rep(1,times=2), 2, rep(1,times=6), rep(3,times=2),
rep(2,times=7)), knots=knots100)## Iteration: 10 (10% completed; 2.29639 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.2222222
## 2 r/delta (overall) 0.7777778
## 3 r/delta (move 1) 1.0000000
## 4 r/delta (move 2) 0.6000000
## 5 r/delta (move 3) NaN## Iteration: 20 (20% completed; 4.26684 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.1578947
## 2 r/delta (overall) 0.4736842
## 3 r/delta (move 1) 0.5555556
## 4 r/delta (move 2) 0.3750000
## 5 r/delta (move 3) 0.5000000## Iteration: 30 (30% completed; 6.21747 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.1034483
## 2 r/delta (overall) 0.4137931
## 3 r/delta (move 1) 0.4166667
## 4 r/delta (move 2) 0.4166667
## 5 r/delta (move 3) 0.4000000## Iteration: 40 (40% completed; 8.30155 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.07692308
## 2 r/delta (overall) 0.46153846
## 3 r/delta (move 1) 0.42857143
## 4 r/delta (move 2) 0.47368421
## 5 r/delta (move 3) 0.50000000## Iteration: 50 (50% completed; 10.31038 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.06122449
## 2 r/delta (overall) 0.38775510
## 3 r/delta (move 1) 0.46666667
## 4 r/delta (move 2) 0.36000000
## 5 r/delta (move 3) 0.33333333## Iteration: 60 (60% completed; 12.33724 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.05084746
## 2 r/delta (overall) 0.40677966
## 3 r/delta (move 1) 0.43750000
## 4 r/delta (move 2) 0.39393939
## 5 r/delta (move 3) 0.40000000## Iteration: 70 (70% completed; 14.61572 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.1304348
## 2 r/delta (overall) 0.4057971
## 3 r/delta (move 1) 0.4500000
## 4 r/delta (move 2) 0.4166667
## 5 r/delta (move 3) 0.3076923## Iteration: 80 (80% completed; 16.61848 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.1265823
## 2 r/delta (overall) 0.3924051
## 3 r/delta (move 1) 0.4545455
## 4 r/delta (move 2) 0.4047619
## 5 r/delta (move 3) 0.2666667## Iteration: 90 (90% completed; 18.54577 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.1460674
## 2 r/delta (overall) 0.3820225
## 3 r/delta (move 1) 0.4000000
## 4 r/delta (move 2) 0.4000000
## 5 r/delta (move 3) 0.3157895## Iteration: 100 (100% completed; 20.49031 secs elapsed)## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.1414141
## 2 r/delta (overall) 0.3535354
## 3 r/delta (move 1) 0.3571429
## 4 r/delta (move 2) 0.3673469
## 5 r/delta (move 3) 0.3181818## The following statement saved the model fits using 100,000 MCMC samples. We won't save
## our fit based on 100 samples. Rather we will load in the results from 100,000.
#save(fit_gvs_knots100,file="bkmr_NHANES_gvs_knots100.RData")
load("./Supervised/BKMR/saved_model/bkmr_NHANES_gvs_knots100.RData")
summary(fit_gvs_knots100)## Fitted object of class 'bkmrfit'
## Iterations: 1e+05
## Outcome family: gaussian
## Model fit on: 2018-12-11 19:32:59
## Running time: 6.46064 hours
##
## Acceptance rates for Metropolis-Hastings algorithm:
## param rate
## 1 lambda 0.08499085
## 2 r/delta (overall) 0.24789248
## 3 r/delta (move 1) 0.33862879
## 4 r/delta (move 2) 0.24140731
## 5 r/delta (move 3) 0.16400060
##
## Parameter estimates (based on iterations 50001-100000):
## param mean sd q_2.5 q_97.5
## 1 beta1 -0.53023 0.04356 -0.61920 -0.44620
## 2 beta2 -0.01427 0.03218 -0.07732 0.04899
## 3 beta3 -0.15517 0.06135 -0.27624 -0.03525
## 4 beta4 -0.01796 0.07069 -0.15548 0.12091
## 5 beta5 -0.06813 0.07909 -0.22397 0.08561
## 6 beta6 -0.08602 0.07921 -0.24138 0.06977
## 7 beta7 0.06920 0.07919 -0.08595 0.22516
## 8 beta8 -0.01479 0.08959 -0.19039 0.16028
## 9 beta9 0.17799 0.10938 -0.03762 0.39132
## 10 beta10 0.04093 0.08263 -0.12039 0.20275
## 11 beta11 0.19093 0.08279 0.02932 0.35303
## 12 beta12 -0.05792 0.03329 -0.12330 0.00741
## 13 beta13 -1.10570 3.06012 -7.08628 4.92819
## 14 beta14 -0.31742 0.74194 -1.76554 1.14300
## 15 beta15 -1.18773 3.39199 -7.81900 5.51197
## 16 beta16 -0.27092 0.78130 -1.80130 1.26774
## 17 beta17 -0.02009 0.16103 -0.33394 0.29761
## 18 beta18 0.05246 0.03224 -0.01091 0.11560
## 19 sigsq.eps 0.74942 0.03396 0.68556 0.81925
## 20 r1 0.00087 0.00688 0.00000 0.01275
## 21 r2 0.00136 0.02158 0.00000 0.01467
## 22 r3 0.00124 0.00812 0.00000 0.01664
## 23 r4 0.00109 0.00624 0.00000 0.01575
## 24 r5 0.00131 0.00832 0.00000 0.01555
## 25 r6 0.00202 0.01103 0.00000 0.02380
## 26 r7 0.00153 0.00783 0.00000 0.01966
## 27 r8 0.00101 0.00778 0.00000 0.01333
## 28 r9 0.00095 0.00613 0.00000 0.01454
## 29 r10 0.00652 0.02102 0.00000 0.05517
## 30 r11 0.00783 0.01333 0.00000 0.04042
## 31 r12 0.00009 0.00166 0.00000 0.00000
## 32 r13 0.00021 0.00481 0.00000 0.00000
## 33 r14 0.00015 0.00230 0.00000 0.00000
## 34 r15 0.01832 0.02523 0.00000 0.07990
## 35 r16 0.00018 0.00230 0.00000 0.00000
## 36 r17 0.00028 0.00350 0.00000 0.00000
## 37 r18 0.00046 0.00590 0.00000 0.00000
## 38 lambda 2.14163 3.51586 0.14823 11.74103
##
## Posterior inclusion probabilities:
## variable group groupPIP condPIP
## 1 PCB74 1 0.41952 0.09320
## 2 PCB99 1 0.41952 0.12710
## 3 PCB118 2 0.86570 0.05656
## 4 PCB138 1 0.41952 0.11709
## 5 PCB153 1 0.41952 0.13816
## 6 PCB170 1 0.41952 0.17263
## 7 PCB180 1 0.41952 0.15098
## 8 PCB187 1 0.41952 0.09926
## 9 PCB194 1 0.41952 0.10159
## 10 PCB169 3 0.61998 0.35298
## 11 PCB126 3 0.61998 0.64702
## 12 Dioxin1 2 0.86570 0.00534
## 13 Dioxin2 2 0.86570 0.00880
## 14 Dioxin3 2 0.86570 0.00737
## 15 Furan1 2 0.86570 0.88037
## 16 Furan2 2 0.86570 0.01047
## 17 Furan3 2 0.86570 0.01243
## 18 Furan4 2 0.86570 0.01867Posterior Inclusion Probabilities (PIPs)
## obtain posterior inclusion probabilities (PIPs)
ExtractPIPs(fit_gvs_knots100)## variable group groupPIP condPIP
## 1 PCB74 1 0.41952 0.093201754
## 2 PCB99 1 0.41952 0.127097635
## 3 PCB118 2 0.86570 0.056555389
## 4 PCB138 1 0.41952 0.117086194
## 5 PCB153 1 0.41952 0.138157895
## 6 PCB170 1 0.41952 0.172625858
## 7 PCB180 1 0.41952 0.150982075
## 8 PCB187 1 0.41952 0.099256293
## 9 PCB194 1 0.41952 0.101592296
## 10 PCB169 3 0.61998 0.352979128
## 11 PCB126 3 0.61998 0.647020872
## 12 Dioxin1 2 0.86570 0.005336722
## 13 Dioxin2 2 0.86570 0.008802125
## 14 Dioxin3 2 0.86570 0.007369759
## 15 Furan1 2 0.86570 0.880374264
## 16 Furan2 2 0.86570 0.010465519
## 17 Furan3 2 0.86570 0.012429248
## 18 Furan4 2 0.86570 0.018666975Data Visualization
Correlation Matrix
### correlation matrix
cor.Z = cor(lnmixture_z, use = "complete.obs")
#pdf(file= "./Supervised/BKMR/saved_model/cor_nhanes.pdf", width=12, height=12)
corrplot.mixed(cor.Z, upper = "ellipse", lower.col = "black")
#dev.off()Model
Trace Plots
### change this for each model you fit and then rerun the code from here to the bottom
modeltoplot <- fit_gvs_knots100 ## name of model object
modeltoplot.name <- "fit_gvs_knots100" ## name of model for saving purposes
plot.name <- "gvs_knots100" ## part that changed in plot name
Z <- lnmixture_z ## Z matrix to match what was used in model
### values to keep after burnin/thin
sel<-seq(50001,100000,by=50)
### access convergence with traceplots
TracePlot(fit = modeltoplot, par = "beta", sel=sel)
TracePlot(fit = modeltoplot, par = "sigsq.eps", sel=sel)
par(mfrow=c(2,2))
TracePlot(fit = modeltoplot, par = "r", comp = 1, sel=sel)
TracePlot(fit = modeltoplot, par = "r", comp = 2, sel=sel)
TracePlot(fit = modeltoplot, par = "r", comp = 3, sel=sel)
TracePlot(fit = modeltoplot, par = "r", comp = 4, sel=sel)
par(mfrow=c(1,1))Univariable Exposure-Response Functions
#### create dataframes for ggplot (this takes a little while to run)
## We will not run these 6 commands that generate all the data for the summaries of
## the exposure-response function. Instead we will read in the previously generated results
## just to save a little time.
#pred.resp.univar <- PredictorResponseUnivar(fit = modeltoplot, sel=sel, method="approx")
#pred.resp.bivar <- PredictorResponseBivar(fit = modeltoplot, min.plot.dist = 1, sel=sel, method="approx")
#pred.resp.bivar.levels <- PredictorResponseBivarLevels(pred.resp.df = pred.resp.bivar, Z = Z,
# both_pairs = TRUE, qs = c(0.25, 0.5, 0.75))
#risks.overall <- OverallRiskSummaries(fit = modeltoplot, qs = seq(0.25, 0.75, by = 0.05), q.fixed = 0.5, method = "approx",sel=sel)
#risks.singvar <- SingVarRiskSummaries(fit = modeltoplot, qs.diff = c(0.25, 0.75),
# q.fixed = c(0.25, 0.50, 0.75), method = "approx")
#risks.int <- SingVarIntSummaries(fit = modeltoplot, qs.diff = c(0.25, 0.75), qs.fixed = c(0.25, 0.75))
#Save the data need to plot the summaries of interest generated
#save(pred.resp.univar, pred.resp.bivar, pred.resp.bivar.levels, risks.overall, risks.singvar, risks.int, file=paste0("./Supervised/BKMR/saved_model/", modeltoplot.name,"_plots.RData"))
# Load in the results, which were computed previously and saved using the command just above this
load(paste0("./Supervised/BKMR/saved_model/", modeltoplot.name,"_plots.RData"))pred.resp.univar %>%
mutate(variable = fct_recode(variable, "PCB 74" = "PCB74",
"PCB 99" = "PCB99",
"PCB 118" = "PCB118",
"PCB 138" = "PCB138",
"PCB 153" = "PCB153",
"PCB 170" = "PCB170",
"PCB 180" = "PCB180",
"PCB 187" = "PCB187",
"PCB 194" = "PCB194",
"1,2,3,6,7,8-hxcdd" = "Dioxin1",
"1,2,3,4,6,7,8-hpcdd" = "Dioxin2",
"1,2,3,4,6,7,8,9-ocdd" = "Dioxin3",
"2,3,4,7,8-pncdf" = "Furan1",
"1,2,3,4,7,8-hxcdf" = "Furan2",
"1,2,3,6,7,8-hxcdf" = "Furan3",
"1,2,3,4,6,7,8-hxcdf" = "Furan4",
"PCB 169" = "PCB169",
"PCB 126" = "PCB126")) %>%
ggplot(aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) +
geom_smooth(stat = "identity") + labs(y = "Estimate", x = "Exposure") +
facet_wrap(~ variable) + theme_bw() +
theme(strip.background = element_rect(fill = "white"))
Bivariable Exposure-Response Functions
pred.resp.bivar %>%
mutate(variable1 = fct_recode(variable1, "PCB 74" = "PCB74",
"PCB 99" = "PCB99",
"PCB 118" = "PCB118",
"PCB 138" = "PCB138",
"PCB 153" = "PCB153",
"PCB 170" = "PCB170",
"PCB 180" = "PCB180",
"PCB 187" = "PCB187",
"PCB 194" = "PCB194",
"1,2,3,6,7,8- hxcdd" = "Dioxin1",
"1,2,3,4,6,7,8- hpcdd" = "Dioxin2",
"1,2,3,4,6,7,8,9- ocdd" = "Dioxin3",
"2,3,4,7,8- pncdf" = "Furan1",
"1,2,3,4,7,8- hxcdf" = "Furan2",
"1,2,3,6,7,8- hxcdf" = "Furan3",
"1,2,3,4,6,7,8- hxcdf" = "Furan4",
"PCB 169" = "PCB169",
"PCB 126" = "PCB126"),
variable2 = fct_recode(variable2, "PCB 74" = "PCB74",
"PCB 99" = "PCB99",
"PCB 118" = "PCB118",
"PCB 138" = "PCB138",
"PCB 153" = "PCB153",
"PCB 170" = "PCB170",
"PCB 180" = "PCB180",
"PCB 187" = "PCB187",
"PCB 194" = "PCB194",
"1,2,3,6,7,8- hxcdd" = "Dioxin1",
"1,2,3,4,6,7,8- hpcdd" = "Dioxin2",
"1,2,3,4,6,7,8,9- ocdd" = "Dioxin3",
"2,3,4,7,8- pncdf" = "Furan1",
"1,2,3,4,7,8- hxcdf" = "Furan2",
"1,2,3,6,7,8- hxcdf" = "Furan3",
"1,2,3,4,6,7,8- hxcdf" = "Furan4",
"PCB 169" = "PCB169",
"PCB 126" = "PCB126")) %>%
ggplot(aes(z1, z2, fill = est)) +
geom_raster() +
facet_grid(variable2 ~ variable1, scales = "free", space = "free",
labeller = labeller(variable1 = label_wrap_gen(5),
variable2 = label_wrap_gen(5))) +
scale_fill_gradientn(colours=c("#0000FFFF","#FFFFFFFF","#FF0000FF")) +
xlab("Exposure 1") +
ylab("Exposure 2") + labs(fill = "Estimate") +
ggtitle("h(Exposure 1, Exposure 2)") + theme_bw() +
theme(strip.background = element_rect(fill = "white"),
panel.spacing = unit(0.05, "lines"),
legend.position = "bottom")
pred.resp.bivar.levels %>%
mutate(variable1 = fct_recode(variable1, "PCB 74" = "PCB74",
"PCB 99" = "PCB99",
"PCB 118" = "PCB118",
"PCB 138" = "PCB138",
"PCB 153" = "PCB153",
"PCB 170" = "PCB170",
"PCB 180" = "PCB180",
"PCB 187" = "PCB187",
"PCB 194" = "PCB194",
"1,2,3,6,7,8- hxcdd" = "Dioxin1",
"1,2,3,4,6,7,8- hpcdd" = "Dioxin2",
"1,2,3,4,6,7,8,9- ocdd" = "Dioxin3",
"2,3,4,7,8- pncdf" = "Furan1",
"1,2,3,4,7,8- hxcdf" = "Furan2",
"1,2,3,6,7,8- hxcdf" = "Furan3",
"1,2,3,4,6,7,8- hxcdf" = "Furan4",
"PCB 169" = "PCB169",
"PCB 126" = "PCB126"),
variable2 = fct_recode(variable2, "PCB 74" = "PCB74",
"PCB 99" = "PCB99",
"PCB 118" = "PCB118",
"PCB 138" = "PCB138",
"PCB 153" = "PCB153",
"PCB 170" = "PCB170",
"PCB 180" = "PCB180",
"PCB 187" = "PCB187",
"PCB 194" = "PCB194",
"1,2,3,6,7,8- hxcdd" = "Dioxin1",
"1,2,3,4,6,7,8- hpcdd" = "Dioxin2",
"1,2,3,4,6,7,8,9- ocdd" = "Dioxin3",
"2,3,4,7,8- pncdf" = "Furan1",
"1,2,3,4,7,8- hxcdf" = "Furan2",
"1,2,3,6,7,8- hxcdf" = "Furan3",
"1,2,3,4,6,7,8- hxcdf" = "Furan4",
"PCB 169" = "PCB169",
"PCB 126" = "PCB126")) %>%
ggplot(aes(z1, est)) +
geom_smooth(aes(col = quantile), stat = "identity") +
facet_grid(variable2 ~ variable1, scales = "free", space = "free",
labeller = labeller(variable1 = label_wrap_gen(5),
variable2 = label_wrap_gen(5))) +
ggtitle("h(Exposure 1 | Quantiles of Exposure 2)") +
xlab("Exposure 1") + theme_bw() + labs(col = "Quantile", y = "Estimate") +
theme(strip.background = element_rect(fill = "white"),
panel.spacing = unit(0.05, "lines"),
legend.position = "bottom")## Warning: Removed 2190 rows containing missing values (geom_smooth).
Overall Mixture Effect
ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +
geom_hline(yintercept=00, linetype="dashed", color="gray") +
geom_pointrange() + theme_bw() +
labs(x = "Quantile", y = "Estimate")
Single Variable
risks.singvar %>% mutate(variable = fct_recode(variable, "PCB 74" = "PCB74",
"PCB 99" = "PCB99",
"PCB 118" = "PCB118",
"PCB 138" = "PCB138",
"PCB 153" = "PCB153",
"PCB 170" = "PCB170",
"PCB 180" = "PCB180",
"PCB 187" = "PCB187",
"PCB 194" = "PCB194",
"1,2,3,6,7,8-hxcdd" = "Dioxin1",
"1,2,3,4,6,7,8-hpcdd" = "Dioxin2",
"1,2,3,4,6,7,8,9-ocdd" = "Dioxin3",
"2,3,4,7,8-pncdf" = "Furan1",
"1,2,3,4,7,8-hxcdf" = "Furan2",
"1,2,3,6,7,8-hxcdf" = "Furan3",
"1,2,3,4,6,7,8-hxcdf" = "Furan4",
"PCB 169" = "PCB169",
"PCB 126" = "PCB126")) %>%
ggplot(aes(variable, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd, col = q.fixed)) +
geom_hline(aes(yintercept=0), linetype="dashed", color="gray") +
geom_pointrange(position = position_dodge(width = 0.75)) + coord_flip() +
theme_bw() +
labs(x = "", y = "Estimate", col = "Fixed Quantile")
Single Variable Interaction Terms
risks.int %>% mutate(variable = fct_recode(variable, "PCB 74" = "PCB74",
"PCB 99" = "PCB99",
"PCB 118" = "PCB118",
"PCB 138" = "PCB138",
"PCB 153" = "PCB153",
"PCB 170" = "PCB170",
"PCB 180" = "PCB180",
"PCB 187" = "PCB187",
"PCB 194" = "PCB194",
"1,2,3,6,7,8-hxcdd" = "Dioxin1",
"1,2,3,4,6,7,8-hpcdd" = "Dioxin2",
"1,2,3,4,6,7,8,9-ocdd" = "Dioxin3",
"2,3,4,7,8-pncdf" = "Furan1",
"1,2,3,4,7,8-hxcdf" = "Furan2",
"1,2,3,6,7,8-hxcdf" = "Furan3",
"1,2,3,4,6,7,8-hxcdf" = "Furan4",
"PCB 169" = "PCB169",
"PCB 126" = "PCB126")) %>%
ggplot(aes(variable, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +
geom_pointrange(position = position_dodge(width = 0.75)) +
geom_hline(yintercept = 0, lty = 2, col = "gray") + coord_flip() + theme_bw() +
labs(x = "", y = "Estimate")